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Abstract 

Recently, a new class of second order Runge-Kutta methods for Ito stochastic diffe- 
rential equations with a multidimensional Wiener process was introduced by Rofiler 
[10]. In contrast to second order methods earlier proposed by other authors, this 
class has the advantage that the number of function evaluations depends only lin- 
early on the number of Wiener processes and not quadratically. In this paper, we 
give a full classification of the coefficients of all explicit methods with minimal stage 
number. Based on this classification, we calculate the coefficients of an extension 
with minimized error constant of the well-known RK32 method [2] to the stochastic 
case. For three examples, this method is compared numerically with known order 
two methods and yields very promising results. 
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1 Introduction 



In recent years, the development of numerical methods for the approximation 
of stochastic differential equations (SDEs) has become a field of increasing 
interest, see. e.g [Iff] and references therein. Whereas strong approximation 
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methods are designed to obtain good pathwise solutions [1], weak approxi- 
mation focuses on the expectation of functionals of the solution. Second or- 
der stochastic Runge-Kutta (SRK) methods for the weak approximation of 
SDEs were proposed by Kloeden and Platen [3], Komori [5], Mackevicius and 
Navikas [6j, Tocino and Vigo-Aguiar [13], and the authors [3]l9] . However, 
these methods were not suitable for problems with high numbers m of Wiener 
processes, because for these methods the number of function evaluations per 
step increases quadratically in m. Recently, new classes of SRK methods were 
introduced by Rofiler PHI?] which overcome this problem. In Section [2] we 
present the one of these classes which is suitable for Ito SDEs. The aim of this 
paper is to give a full classification of all the explicit methods within this class 
with minimal stage number, which is done in Section |3j As an application, 
in Section H] we extend the well known RK32 scheme [2] to an SRK method 
with minimized leading local error term. The performance of this method is 
illustrated by some numerical examples in section 

We denote by (X(t)) t ^i the solution of the rf-dimensional Ito SDE defined 



with an m-dimensional Wiener process (W(t)) t >o and / = [io>^1- We assume 
that the Borel-measurable coefficients a : I x R d — > R d and b : I x R d -> 
M dxm satisfy a Lipschitz and a linear growth condition such that the Ex- 
istence and Uniqueness Theorem [3] applies. In the following, let V>(t,x) = 
(b l, i(t,x))i<i<ti E R d denote the jth column of the diffusion matrix b(t,x) for 
j = l,...,m. 

Let a discretization Ih = {to, ti, • • • , t?v} with to < h < . . . < = T 
of the time interval I = [to,T] with step sizes h n = t n+ \ — t n for n = 
0, 1, . . . , AT — 1 be given. Further, define h = max < n <Ar h n as the maximum 
step size. Let C l P (R d , R) denote the space of all g E C l (R d , R) fulfilling a poly- 
nomial growth condition and let g E Cp\l x R d , R) if g(-,x) E C k (I,R) and 
g(t, •) E C l P (^ d , K) for all t E I and x E R d % 



Definition 1.1 A time discrete approximation Y = (Y(t)) te j h converges weakly 
with order p to X as h — )■ at time t E Ih if for each f E C 2 p p+1 ^ (R d , R) exist 
a constant Cf and a finite 5q > such that 



by 



dX(t) = a(t, X(t)) dt + bit, X(t)) dW{t), 



X(to) = %o, 



(1) 



E(f(X(t)))-E(f(Y(t)))\<C f W 



(2) 



holds for each h e]0,5q[. 
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2 Stochastic Runge-Kutta methods 



We consider the stochastic Runge-Kutta methods introduced in [10] for the 
weak approximation of SDE ([1]). Therefore, the (i-dimensional approximation 
process Y with Y n = Y(t n ) of an explicit s-stage SRK method is defined by 
Y = xq and 

s 

Y n+1 = Y n + E a i a (tn + c\°'h n , Hj®) h n 

i=l 
s m 

+ EE^ ) b k (t n + c^h n ,H^)i {k) 

i=l k=l 

s m j 

+ J2J2^ ) b k (t n + c^h n ,Ht ) ) 1 -^ (3) 
i=l k=l v ftn 

s m 

+ Y,Y,ti 3) b k (t n + c? ) h n ,H? ) )i ik) 

i=l fc=l 

s m 

i=i /c=i 

for n = 0, 1, . . . , N — 1 with stage values 

tff = Y n + E Ag } o(t» + ef #f ) fr n 
i— 1 m 

i=i J=l 

3=1 

+ E^; ) 6 fc (t n + 4 1) / in ,i/f ) )vt 

i=l 

tff 3 = F n + E 4? + c f ^, ^f) 

i=l 

s m f 

+ EE4 ) ^ + 4 1) ^X ) )%^ 

j=l 1=1 V^n 
l^k 

for i = 1, . . . , s and jfe = 1, . . . , m. Here, a, . . . , /3 (4) , c w G M s and A®, 
B ( q ) e Rsxs f or o < g < 2 with A§ } = £§ } = for j > i and < q < 1 
are the vectors and matrices of coefficients of the SRK method, c^ q > = A^e 
for < q < 2 with a vector e = (1, . . . , 1) T . In the following, the product 
of column vectors is defined component-wise. The coefficients of the SRK 
method ([3]) are determined by the following Butcher tableau: 
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c (0) 




5(0) 




c (l) 




5(1) 




c (2) 


^(2) 


£(2) 






T 

a 










/3( 3 ) r 


/3W T 



are three-point distributed random variables with P(I(k) = ±V3 h n ) = | 
and P(/(fc) = 0) = |. Further, 1^ are defined by 

{l(i(k)i(i) - VK,i(k)) ftk<i 
\{i(k)I(i) + VKI(i)) iil<k (4) 
Wl)-h n ) Hk = l 

with two point distributed random variables 1^) satisfying P(I(k) — iv^n) = 

1 1 

2 ' 

By the application of the multi-colored rooted tree analysis [8], order con- 
ditions for the coefficients of the SRK method ([3]) can be easily determined. 
As a result of this, the following Theorem 12.11 due to Rofiler [10] gives order 
conditions for the SRK method ([3]) up to order two. 

Theorem 2.1 Let a\b iJ G C 2 /(I x R d , R) for 1 < i < d, 1 < j < m. If the 
coefficients of the SRK method (TJ|) fulfill the equations 

1. « T e = l 2. /?( 4 ) T e = 3. /^e = 

4. (/3« T e) 2 = l 5. (3^ T e = 6. ^W T S«e = 

7. /3( 4 ) T A( 2 )e = 8. pW T BMe = 9. /^(E^e) 2 = 

then the method attains order 1 in the weak sense. In addition, if a 1 , e 
Cp' 6 (J x R d , R) for I <i < d, 1 < j < m and if the equations 



10. 


a T A^e = | 




11. 


a T {B^ef = \ 




12. 


(^ T e)(a T B^e) = 


1 

2 


13. 


(^ T e)(^ T A^e) 


i 

2 


14. 


pm T A m e = o 




15. 


pm T B m e = i 




16. 


/3 w T jB (2) e = i 




17. 


{^ T e)(^ T (B^e 




18. 






19. 


/g (i) T ( jB (i) (jB (i) e) ) = 





20. 


^(3) T (jB (2) (jB (l) e)) = 





21. 


(3^ T (B^(B^\B^ 


e)))=0 


22. 


^(D T (A (i) (jB (o) e)) = 





23. 


^(3) T (A (2) (jB (o) e)) = 
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94 


/ o(4) T M (2) x2 _ q 






/Q(4) T (' / l(2)/ / i(0) (3 ^ _ n 

y(J 1/1 1 il G 1 1 U 


9fi 


n T '( /?(°) r - n 

Cl I J_> I XJ t- 1 1 — u 




97 


/9(2) T >i(i)p _ n 

/J i 1 o — U 


98 


^ (1)T rf/4 (1) pVR (1) pH — n 






/3(3) T // 4(2) Wn(2) \N _ n 


ou. 


B^ T ( A^(E^p)) - o 

fj 1/1 I XJ G 1 1 U 




O 1 . 


fj \ i\ \ XJ C J J — u 


^9 


p \\ n e ) l^ 1 e )) — u 




OO. 


/ o(4) T M (2)/ R (0) n2n _ n 




/ o(2) T M (l)/ R (0) x2n _ n 




oo. 


/Q(1) T |'r(1)M( 1 V^ - n 


OU . 


o(3) T ( B (2) ( a(1) \\ _ q 




O I . 


B (2) T ( B (1) \2 _ q 
/X I XJ O J — u 


38 
oo. 


R^) T (B {2) (B {1) eW - 

|X I XJ I XJ o 1 1 — \J 




39 


B i2)T (B w (B {1) p)) - 


in 


RW T (n^p] 3 - o 




41 


/ o(3) T rR (2) x3 _ n 




yLV 1 XJ I U CI 1 \J 




43 
'-to. 


fl(3) T /n(2) ( ' jB (l) N2X _ q 




ft(4) T / R (2) \4 _ r, 

/J 1 XJ c- J — U 




4^ 


/ o(4) T / R (2)/ R (l) xn2 _ r, 


46 


8^ T ((B^p)(B^(B^p))) - 
P \\ n e )\ n \ n e ))) — 


n 

u 


47 


ry T ((B^p)(B ( -°UB^p))) - Q 




P U 71 \ n e ))\ n e >) — 


n 

u 


4Q 


o(3) T / M (2)/t ? (0) xw R (2) _ n 


50 


R^ T (A^(R^(R^p))) - 
P I 71 1-° 1-° v))) ~ u 




^1 
tj X . 


o(3)T ( *(2)/d(0)^d(1) \\\ _ Q 
/j i^/i i^iD i^.d e;;; — u 


£ o 

oz. 


p v y U-D (Jl 'e))\n y 'e)J = 


n 
U 


Oo. 




54. 


/3 (3) T (S (2) (A (l) (jB (0) e))) =Q 




55. 


(3 {1)T ((B^e)(B^(B^e))) = 


56. 


/3( 3 ) T (( J B( 2 ) e )( J B^( J B (1) e))) = 





57. 


/g (i) T (jB (i) (jB (i) (jB (i) e))) = 


58. 


^ T ((B^e)(B^(B^e) 2 ))-- 


= 


59. 


/3( 4 ) T (( J B( 2 )e)( J B( 2 )( J B«( J B( 1 )e)))) 



are fulfilled, then the SRK method |3j] attains order 2 in the weak sense. 



It turns out that explicit order one SRK methods need at least s = 1 stage 
while order two SRK methods need s > 3 stages. This is due to e.g. conditions 
4., 6. and 17., which can not be fulfilled in the case of s < 2 stages for 
explicit order two SRK methods. In the following, we distinguish between the 
stochastic and the deterministic order of convergence. Let Ps = P denote the 
order of convergence of the SRK method if it is applied to an SDE and let 
Pd with pb > ps denote the order of convergence of the SRK method if it is 
applied to a deterministic ordinary differential equation (ODE), i.e., SDE ([1]) 
with 6 = 0. We also write (pd,Ps) m the following. 
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3 Parameter families for SRK methods 



3. 1 Coefficients for SRK methods of order (1,1) 



First, we analyze explicit SRK methods fl3]) of order po = p$ = 1 with s — 1 
stage. Considering the order one conditions 1.-9. in Theorem 12. 1[ the corre- 
sponding coefficients are uniquely determined for ci G { — 1,1} by 

«i = l, & (1) = ci, (3{ 2) = 0, (3[ 3) = 0, M 4) = 0. (5) 

The resulting class of SRK schemes coincides with the well-known Euler- 
Maruyama scheme. 



3.2 Coefficients for SRK methods of order (2,1) 



Next, we consider the case of s = 2 stage explicit SRK methods (j3J). As already 
mentioned in Section [2], it is not possible to attain order ps = 2. However, we 
can find some SRK methods of order pp = 2 and ps — 1 corresponding 
to the following parameter family: From condition 1. of Theorem 12.11 follows 
«i = 1 — «2 and taking into account the order two condition 10. we obtain 
a 2 = ^T(oT f° r ^21 7^ 0. Further, condition 2. yields = —(3^, condition 

J (3) _ o(3) 



2,1 



3. results in f3{ = —f3 2 
condition 4. holds for 
condition 6. we need that 



}( 2 ) 



-4 2) while 



j (3) 



and condition 5. is fulfilled if f3\ 

= Ci 

0, considering condition 8 



02 with Ci G {—1,1}. Finally, considering 



or B^ + B^ 



or B$ 



(2) (2) ■ — 

B21 + B22 and for condition 7. and 



analogously that /3 2 

9. that (3 { 2 ] = or 4i + 4? = 4i +4? and + 5g } ) 2 
hold. Thus, this class of SRK methods is determined by 



/ R (2) 
(#21 



#22 , 



a: 

A® 
B (o) 



1 - 



1 



1 



2c 2 2c 2 





c 2 

0' 

c 3 



c 16. 



ci t |_— 1, ij- ana c 2 , . . . , c 17 
= 0, c 7 (c 9 + cio - cn - C12) 



Ci — C4 C4 

"0 

c 8 

" 

Cl7 



pw T = 


C5 -C5 




c 7 -c 7 


A® = 


Cg Cio 




Cll C12 


B^ = 


C13 C14 




C15 Ci6 


7 = 0, C 6 (ci 3 + Ci4 


f ~ (C15 


+■ ci 6 ) 2 ) 



(6) 
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3.3 Coefficients for SRK methods of order (2,2) 



Now, we consider explicit SRK methods of order pa = ps = 2 with s = 3 
stages. Then, the SRK schemes of the class under consideration are completely 
characterized by the following families of coefficients which follow from the 
order conditions in Theorem l2.lt 

We have ai = 1 — «2 — ot 3 due to condition 1. From condition 4. it follows that 
fiP = ci - /3 { 2 1] - I3[ 1] with ci G {-1, 1}. Due to conditions 2., 3., 7., 24., 16., 
14., 18. and 8. we need £f=i 4? = £- =1 4? = Ef=i4?- From conditions 
3., 8., 18. and 41. follows that £f J=1 -Bff = and that k := Y% =1 B$, 

(3) (3) 

i = 1, 2, 3 are pairwise different. Further, we have f}\ = 2(b 1 -b 2 )\b 2 +2b 1 ) ' <^2 = 
2(6 2 -6 1 )(6 1 +26 2 ) > ^ = 2{2b,+bl)(b 1+ 2b 2 ) • From conditions 2., 9. and 16. we obtain 

o(4) _ &i o(4) _ 6 2 o(4) _ -61-62 WUh 44 if 

1 1 (61 -62) (62+261)' ^2 - (62 -61) (61 +26 2 )' ^3 (26i+6 2 )(6i+26 2 )- VV L 1 

follows now from the above that exactly for one i from 1, 2, 3 it holds bi = 0. 
Without loss of generality we can assume that b\ = 0. Due to conditions 15. 
and 37. we need B$ ^ and from conditions 17. and 40. follows $3 ^ 0. 
Now, by condition 19. follows that B^ 2 = and we deduce from 6., 17. and 40. 
that B$ = -B { 2i- With conditions 5., 15. and 37. follows that f3 { 2 2) = ^y, 

= -^tt and (3{ 2) = 0. Conditions 4., 6. and 17. yield =c x - 

and P 2 = 03 = 4(^(1) ) 2 • F rom condition 56. and 58. we obtain B^ = —B^ 

(2) (2) 21 (2) (2) 

and B 23 = —B 33 . Condition 43. yields now B y y> = -B\l'. Condition 46. gives 

(2) (2) (2) 

B 32 = B 33 . Condition 20. leads to B 13 = 0. Now, due to 13. and 28. we need 
that 4¥ = {B^f and A« = (£«) 2 - A 32 . 

To fulfill 11., 12., 22., 23., 30. and 33., we obtain the following cases: 

1) Bg> = 0, + = ci, a 3 = |, 4? = 4? = 4?, 

2) A« = 0, B$? = ci, B$> + Bg> = 0, a 2 = |, 4? = 4? = 4?, 

3) 42 } = 0, Bg> = 2$ + 2$ = c 1; a 2 + « 3 = |, 4? + 4? = 4? + 4? = 
4(2) , 4(2) 

■ / *32 ' -^33 5 

4 A 4W _ n o(0) / q / o(0) „(0) / d(0) 4 (2) _ ,(2) ,(2) _ ,(2) .(2) _ 
^32 — u > -°21 7^ u T D 3\ + -°32 7^ -°21 5 ^22 ~ ^32 > ^23 _ ^33 > A \2 ~ 

A (2) ,( A (2) A (2)^ B^+B^ i l-gjB^+B^) i 

^32+1^33 ^13 J R (o) ' a 2 — 2 R (0)/ R (0) R (0) „(o)o a 3 — 2 , (o) R (0) VR (0) R (0) R (0) 

^21 -^21 Wzi ""^31 _±f 32 J 1^31 +^32 A^l _±f 31 _±f 32 

However, from 26. and 51. it follows 

a) b£> = or 

b) «3 = 0,4 2 3 ) + 4 2 3 ) = 24 2 3 ) . 

Finally, the equations 10. and 25. imply the cases 
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i) a*{Ag>-A)jg)?n.. i (A%-A 

| (2) 



|(2)> 



I (2) 



A^'-A (2) 

^22 ^32 



■2 Q2 (^-^)-<*3(^-^) 



(2). 
32 , 

A 



A (P) 



^23 ^33 



2a 2 (A$-A$)-a 3 (A$ 



A 



(0) 



(0) 



ii) 4? 



*23 



33 5 ^22 



32 ; 

22 ^32 

(2) „. / n ,(0) 



i-2 Q 3( 4°>+4°>) 



m 



.4 



(2) 
23 



A (2) .(2) 
^33 ) ^22 



4 2 2 } ,« 2 = 0, « 3 ^0, 4l } = 2^ 



A 



(0) 
32 • 



With these settings, all the remaining order conditions are now fulfilled. 
Summarizing our results, we have the following classification for the SRK 
schemes of order pr> = ps = 2 for the considered class with s = 3 stages: For 
Ci G { — 1,1} and c 2 , c 3 , c 4 , c 5 G R with c 3 7^ and c 4 7^ holds 



/3( 4 ) T 

B (2) 







ci 
2c 2 



ci 



-3 ™-3 ^"3 

1 1_ 

2c 4 2c 4 

0' 


C3 - c 2 c 2 



C4 + 2C5 — c. 
-c 4 - 2c 5 c 5 





-c 5 

c 5 



A. 

2c 3 



1 

"2c 3 



2d 



ci 



ci 



0' 

c 3 
-c 3 



(7) 
(8) 

(9) 
(10) 



Now, the following cases are possible: 

In the case [W| we get with eg, . . . , ci 2 el that 



a 



B (o) 






ci 



(11) 



A (o) 



0' 


en 1 - C12 



A (2) 



C 6 

c 6 
c 6 ■ 



- c 7 c 7 c 8 

" Cg C9 Cg 
ClO C10 Cg 



In the case [lMi|) we obtain with c 6 , . . . , c i2 G R and Ci 7^ that 



(12) 



\ — C10 C10 I 



£(0) 





d 



(13) 



A (0) 





l-cii 
2cio 

Cll - Ci 2 C12 



0' 




A (2) 



c 6 — C 7 C7 Cg 
C 6 - Cg Cg Cg 
C 6 - Cg Cg Cg 
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Considering the case nMIT|) we obtain with cq, . . . , en G M that 



a 



1 i 

2 u 2 



For the case 
coefficients 



"0 

A® = cio 
1 - c a en 

we get with C2 = 



B (o) 



A (2) 



in 





ci 



c 6 - c 7 c 7 c 8 

c 6 ~~ c 9 c 9 c 8 

c 6 ~~ c 9 c 9 c 8. 

and c 6 , . . . , C12 G ' 



1 „ i 

2 _ Cl1 2 



Cll 



B (0) 



In the case 



"0 0' 

A (0) =10 
Ci2 -cia 
we obtain with C2 



C 6 

c 6 • 




ci 


- c 8 c 7 c 8 

- c 9 c 7 c 9 

ClO C 7 Cxo 



(15) 
(16) 

c 9 7^ cio the 
(17) 

(18) 



in © and c 6 , . . . , Ci 2 6 K the coefficients 



A (o) 



For the case 
coefficients 



' 1 ' 1 \ — Cio | Cio 



0~ 

1 - 2ci cii 

Cll - C12 C12 0_ 

we get with C2 = in 



B (o) 




ci 




A (2) 



c 6 ~ c 8 c 7 c 8 

c 6 - c 9 c 7 c 9 
c 6 - c 9 c 7 c 9 

and Cq, . . . ,Cn £ I 



1 i 

2 2 u 



S (0) 



In the cases 

C6, . . . , C12 G 




1 

Cll 

a 



0" 

, A^ = 

en 0_ 

I]) respectively 





Cl 
Cio 

C6 - C 8 
C6 - C 9 





-c 10 

C7 

c 7 



c 8 
Cg 



(19) 
(20) 

C 8 7^ c 9 the 
(21) 

(22) 



c 6 ~ 2C 8 + Cg C 7 2c 8 — Cg_ 

T]) we obtain with c 2 = in (EI) and 



o 



l l 



5(0) 




ci 

Cg - Cio Cio 



A (0) 



0' 

1 

Cll C12 



c 6 c 7 c 8 
c 6 c 7 c 8 
c 6 c 7 c 8 



(23) 



(24) 
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where eg = in the case , eg = c\ for I3UHI) and 7^ eg 7^ c\ for I41dtij) , 
respectively. 

Considering the case I3fefi| we find with eg = in and c 6 , . . . , Ci 2 G R, 
c 9 7^ c 10 that 

"0 01 



a 



1 1 

2 2 



Cll Cn 



S (0) 



A (0) 



'0 00 
1 

1 - C12 C12 o_ 
For the case EMI)) we obtain with c 2 = in 
that 



ci 

ci 

C6 C7 C 8 

C 6 Cg C 7 + C 8 - Cg 

C6 C10 C 7 + C 8 - C10 

HD and c 6 , . . . , C12 G 



(25) 
(26) 

c 10 7^ 2' 



















"0 


0" 








T 

a = 


1 

2 


1 

2 


ClO 


ClO 






Cl 









(27) 





































0" 








"c 6 


C7 




c 8 






1— 2ciocn 












c 6 


cg 


c 7 + c 8 


- Cg 


• (28) 


l-2cio 




1 


.Cll — C12 


C12 


0_ 








c 6 


eg 


c 7 + c 8 


- Cg_ 





A (0) 



In the case lSSnj) respectively ffljn]) we get with C2 = in (jUJ) and C6, 
with 7^ C10 



C12 g 



A (o) 



0' 

cn 
1 - C12 C12 



A m 







" 








I 




ClO 












_ c l 








"c 6 + (eg - c 7 )(l 







(29) 



c 8 + (c 9 -c 7 )^ C- 



<-'■(-> 

C 6 



C 8 

c 8 



eg 
eg 
(30) 



where cio = ci in the case l3Mhj) and cio 7^ ci in the case l4mnj) . 
Considering the case l3fcHl) we obtain with C2 = in ([9]) and eg, . . . , c\\ G 
c 7 7^ c 9 , that 



T 



a 



1 i 

2 2 u 



A (0) 



0' 

1 

1 - cn cn 



ci 

Cl — Cio Cio 0. 

c 6 c 7 c 8 

C 6 Cg C 7 + C 8 - Cg 

C 6 2c 7 - Cg Cg + C 8 - C 7 



Next, we have the case Batil) with C2 = in (jUJ) and C6, . . . , C13 G 

7^ cio 7^ cn 7^ 0, cn 7^ ci, 





T 

a 




1 - 


1— ei(eio+cii' 
2ci cn 




' 





0" 






ClO 








, A^ = 




Cll 











2 cio(cio-cn) 



1 1-cicin 

2 cn (cio— cn) 



0' 

cio cil(cil— cio)— ci2(l— cicip) g g 
cn cicii-1 

Cl2 - C13 C13 0. 



(31) 

(32) 
with 

(33a) 
(33b) 



Published in Applied Numerical Mathematics 59 (2009) no. 3-4, pp. 582-594, \db~i: 10.1016/j.apnum.2008.03.012 



C 6 + (c 9 



c 7 )(l-2n) 



.„, C 8 + (Cg-C 7 )^ C 7 

A< 2 > = C 6 C 8 Cg . (33c) 

C6 Cg Cg_ 

In the remaining cases IHty -niHEi i , IMnj) , IMn|) JMn]) , MJ) , HE]) and HE]) there 
doesn't exist a solution. 



4 Application: An SRK scheme with minimized error coefficients 



Based on the classification given in sectional as an example we will now extend 
the well known method RK32 of Kutta [2] to an SRK method of order (3,2). 
The Butcher array of RK32 is obtained from family ( 13"3"|) by setting c w = 6= ^q , 

c n = 3± 5 ) c i2 = 1, C13 = 2. Due to some symmetry in the method, the sign 
of Ci has no influence, and we choose c\ = 1. Now, we want to determine 
the remaining coefficients by minimizing the expectation of the local error. 
Therefore, we distinguish between the cases m = 1 (only one Wiener process) 
and m > 1. In the case of m = 1, to save computational effort we require that 
equals the zero matrix, because then we don't have to evaluate the stages 
#(*). 

Now, let lef(h) be the weak local error of the method starting at the point 
(t, x) with respect to the functional / and step size h, i. e. 

le f (h) = E (f(Y(t + h)) - f(X(t + h))\Y(t) = X(t) = x). 

As in the deterministic case, by the colored rooted tree analysis one obtains 
the representation 

le f (h) = lec t F ( t )( x ) h3 + °( h4 )> 

teTS(A) 

p(t)=3 

where TS(A) denotes a set of trees, p(t) the order of the tree t, F(t) the ele- 
mentary differential connected with the tree t and lec t a coefficient depending 
only on t and the numerical method (see [8. 9|T0] for details). 
Let lec = (/ect)teTS(A) be the vector of these coefficients. In the following, 
we want to minimize ||Zec|| in the Euclidean norm. Then, using again the 
rooted tree analysis, a tedious calculation (for m = 1 there exist 164 rooted 
trees of order three) yields that in the Euclidean norm we have in the case of 

_ 6-V6 
c io - —LO- 
ui i|2 _ 60500644673+24530366872 y / 6-(217+88v / 6)(128250000c§-92062500c|) 
'I eC H _ 24000000(24+11^6)2 

~38 



which is minimized by C3 = ±3y^-, which gives ||/ec|| ~ 1.275. In the case of 

c io — 6+ io® instead, we would obtain ||/ec|| ~ 1.296, so we choose the minus 
sign in the following. The remaining coefficients of the method, C4 and C5, 
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-1 2 



e-Ve 
10 

3+2^ 
5 



:S 1.2 
491 

342 
491 




214 
513 

214 
513 



1105 
991 



1105 
991 



491 
513 

491 

513 



221 
4955 



221 
1955 



491 
513 

491 
513 



221 
4955 



221 
4955 



193 
(isl 



491 
1368 



491 
1368 



O P. \ / Q O 




4955 
7072 



4955 
14144 



Table 1 

Coefficients of the SRK scheme DRI1 with p£> 



4955 
14144 



3 and ps = 2. 



are determined by considering ||Zec|| in the case of two Wiener processes. Its 

minimal value 2.859 is attained for c 5 = =tf?|\/j|^ c 4 — -pW J 2 ^- For c 3 , c 4 
and C5, the method is invariant to the choice of the sign, so we obtain finally 
the scheme DRI1 presented in Table [U 

In the case m > 1, we cannot avoid completely the evaluation of the stages 
by letting equal to the zero matrix, so one could try to use the additional 
degrees of freedom to minimize the local error. The resulting method differs 
from DRI1 only in A^ 2 \ which is now given by 



' 2(-3442595658+1259007085y / 6) 
1554073317(-6+v / 6) 

!(7-2^6) 
£(7-2>/6) 



8(212963260+73915807^6) 
1554073317(-6v / 6) 



(4(-1111473969+371403611 



23311099755 



81 
_8_ 

81 



(3 + V6) 



81 
J_ 

81 



-3 
-3 



V6) 



and by c*- 2 -* which fulfills now c*- 2 -* = (|, |, |). For m—1, this method has again 
1 1 Zee 1 1 ~ 1.275. But in the case of m = 2, we achieve ||/ec|| ~ 2.765. However, 
this is only 3.4% less than ||Zec|| achieved by DRI1. Due to the additional 
m function evaluations needed for 7^ in the case m > 1 (because for 
A^ = we would have = H[ k \ k = 1, . . . , m), we favor the SRK method 
DRI1 also for m > 1. 
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5 Numerical example 



In the following, the SRK scheme DRI1 presented in Section H] is applied to 
three test equations in order to analyze its order of convergence in compar- 
ison to some well known schemes. Therefore, the functional u = E(f(X(t))) 
is approximated by a Monte Carlo simulation. The performance of DRI1 is 
compared to the second order SRK schemes PL1WM due to Platen pQ, NON 
due to Komori [5], which in contrast to all other schemes is designed for the 
weak approximation of Stratonovich SDEs, RDI3WM and RDI4WM due to 
the authors [3] and the extrapolated Euler-Maruyama scheme EXEM [12] also 
attaining order two, which is given by 2 E(/(Z /l//2 (t))) — ~E(f(Z h (t))) based on 
the Euler-Maruyama approximations Z h l 2 {t) and Z h {t) calculated with step 
sizes h and h/2. The sample average UM,h = jj Hk=\ fO^i^i u k)), ^>k G of 
M independent simulated realizations of the considered approximation Y(t) 
is calculated in order to estimate the expectation. In the following, we denote 
by [L — um,h — E(/(X(£))) the mean error and by a 2 the empirical variance of 
the mean error. Further, we calculate the confidence interval with boundaries 
a and b to the level of 90% for the estimated error fi (see [lj for details). 

As first example, we consider the non-linear SDE [4"f6] 

dX(t) = (^X(t) + ^X(t) 2 + dt+y/X(t) 2 + ldW(t), X(0) = 0, (34) 

on the time interval I = [0, 2] with the solution X(t) = sinh(t + W(t)). Here, 
we choose f(x) = p(arsinh(x)), where p(z) = z 3 — 6z 2 + 8z is a polynomial. 
Then the expectation of the solution can be calculated as 

E(f(X(t))) =t 3 - 3t 2 + 2t . (35) 

The solution E(f(X(t))) is approximated with step sizes 2 _1 ,...,2 -4 and 
M = 10 9 simulations are performed in order to determine the systematic error 
of the considered schemes at time t = 2. The results for the applied schemes 
are presented in Table [2j The orders of convergence correspond to the slope 
of the regression lines plotted in the left hand side of Figure Q] where we get 
the order 1.80 for EXEM, order 1.81 for PL1WM, order 1.93 for RDI3WM, 
order 2.01 for RDI4WM, order 2.41 for NON (applied to the corresponding 
Stratonovich version of fl34l ) and order 2.01 for the scheme DRI1. 
Of course, these results have to be related with the computational effort of the 
schemes which we take in the following as sum of the number of evaluations of 
the drift function a and of each diffusion function V , 1 < j < m, as well as the 
number of random variables that have to be simulated. Then we can compare 
the computational effort versus the errors of the analyzed schemes. The results 
are presented in the right hand side of Figure [TJ The Platen scheme, RDI4WM 
and the new scheme DRI1 yield comparable results and all three are better 
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Table 2 

Mean errors, empirical variances and confidence intervals for SDE (|34p . 





h 




*2 


a 


b 




2~ 


-i 


-1.359E-00 


2.990E-06 


-1.359E-00 


-1.359E-00 


T?YT?l\/r 
tVA.il/ JY1 


2" 

2" 


-2 
-3 


-6.614E-01 
-1.945E-01 


7.315E-06 
8.629E-06 


-6.620E-01 
-1.952E-01 


-6.607E-01 
-1.938E-01 




2" 


-4 


-5.570E-02 


9.014E-06 


-5.641E-02 


-5.499E-02 




2" 


-1 


-3.837E-01 


1.885E-06 


-3.841E-01 


-3.834E-01 


PT 1 Wfl\/f 
rijl VV 1V1 


2" 

2" 


-2 
-3 


-1.165E-01 
-3.348E-02 


3.207E-06 
2.475E-06 


-1.169E-01 
-3.386E-02 


-1.161E-01 
-3.311E-02 




2" 


-4 


-8.949E-03 


3.447E-06 


-9.390E-03 


-8.509E-03 




2" 


-1 


-3.926E-01 


1.400E-06 


-3.929E-01 


-3.923E-01 


SXLJLo VV 1V1 


2" 
2" 


-2 
-3 


-1.041E-01 
-2.748E-02 


2.787E-06 
2.427E-06 


-1.045E-01 
-2.785E-02 


-1.037E-01 
-2.711E-02 




2" 


-4 


-7.054E-03 


1.813E-06 


-7.373E-03 


-6.734E-03 




2" 


-1 


-3.760E-01 


1.488E-06 


-3.762E-01 


-3.757E-01 


rlU14 VV 1V1 


2" 
2~ 


-2 
-3 


-9.454E-02 
-2.318E-02 


2.823E-06 
2.441E-06 


-9.494E-02 
-2.355E-02 


-9.414E-02 
-2.281E-02 




T 


-4 


-5.816E-03 


1.816E-06 


-6.135E-03 


-5.496E-03 




T 


-1 


-3.393E-01 


2.530E-06 


-3.396E-01 


-3.389E-01 


NON 


2~ 


-2 


-4.354E-02 


3.371E-06 


-4.398E-02 


-4.311E-02 


2" 


-3 


-9.707E-03 


2.208E-06 


-1.006E-02 


-9.355E-03 




2" 


-4 


-2.119E-03 


2.952E-06 


-2.526E-03 


-1.711E-03 




2" 


-1 


-3.684E-01 


1.720E-06 


-3.687E-01 


-3.681E-01 


DRI1 


2" 
2" 


-2 
-3 


-9.271E-02 
-2.270E-02 


2.939E-06 
2.122E-06 


-9.312E-02 
-2.304E-02 


-9.231E-02 
-2.235E-02 




2" 


-4 


-5.617E-03 


2.931E-06 


-6.023E-03 


-5.212E-03 




Id h ld(function and random number evaluations) 



Fig. 1. Orders of convergence and computational effort per simulation path versus 
precision for SDE (f34"|) . 

than RDI3WM and much more efficient than the extrapolated Euler method. 
For higher precision, NON performs best. 

As a second example, a multi-dimensional SDE with initial value X(0) = 
(1, 1) T and noncommutative noise driven by a 2-dimensional Wiener process 
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Table 3 

Mean errors, empirical variances and confidence intervals for SDE (|36p . 





h 


A 


*2 


a 


b 


EXEM 


2-u 

2- 1 
2- 2 
2- 3 


-2.165E-05 
-7.684E-06 
-2.266E-06 
-6.078E-07 


1.678E-14 
1.418E-14 
1.501E-14 
3.567E-14 


-2.169E-05 
-7.724E-06 
-2.308E-06 
-6.725E-07 


-2.160E-05 
-7.643E-06 
-2.224E-06 
-5.432E-07 


PL1WM 


2-u 

2- 1 
2- 2 
2- 3 


3.093E-05 
4.947E-06 
1.071E-06 
2.435E-07 


9.082E-15 
1.085E-14 
5.886E-15 
4.652E-15 


3.090E-05 
4.906E-06 
1.041E-06 
2.172E-07 


3.097E-05 
4.987E-06 
1.101E-06 
2.699E-07 


RDI3WM 


2-u 

2- 1 
2-2 

2- 3 


-1.092E-05 
-2.335E-06 
-5.143E-07 
-1.285E-07 


2.481E-15 
8.234E-15 
5.519E-15 
4.581E-15 


-1.094E-05 
-2.370E-06 
-5.431E-07 
-1.546E-07 


-1.090E-05 
-2.299E-06 
-4.856E-07 
-1.023E-07 


RDI4WM 


2-u 

2- 1 
2-2 
2- 3 


-9.312E-06 
-1.893E-06 
-4.096E-07 
-1.035E-07 


3.403E-15 
8.765E-15 
5.591E-15 
4.597E-15 


-9.334E-06 
-1.929E-06 
-4.386E-07 
-1.297E-07 


-9.289E-06 
-1.857E-06 
-3.807E-07 
-7.724E-08 


NON 


2-u 

2- 1 
2-2 
2- 3 


6.396E-06 
1.548E-06 
3.799E-07 
8.544E-08 


1.588E-14 
1.266E-14 
6.172E-15 
4.713E-15 


6.347E-06 
1.504E-06 
3.495E-07 
5.889E-08 


6.445E-06 
1.591E-06 
4.102E-07 
1.120E-07 


DRI1 


2-0 

2- 1 
2-2 
2- 3 


-9.391E-06 
-1.908E-06 
-4.127E-07 
-1.041E-07 


3.332E-15 
8.710E-15 
5.587E-15 
4.597E-15 


-9.414E-06 
-1.944E-06 
-4.416E-07 
-1.304E-07 


-9.369E-06 
-1.872E-06 
-3.838E-07 
-7.792E-08 



is considered [3]: 

x i \_(-m ° ^ ( xl ) dt+ ( ^ xl ^ xl )d( wl ) 

(36) 

Here, we are interested in the second moments which depend on both, the drift 
and the diffusion function (see [I] for details). Therefore, we choose f(x) = 
(x 1 ) 2 and obtain 

E(f(X(t))) = exp(-t) . (37) 

We approximate E(f(X(t))) at t = 10 by M = 8 ■ 10 7 simulated trajectories 
with step sizes 2°, . . . ,2 -3 . The results for the schemes in consideration are 
presented in Tableland Figure [2J Here, the order of convergence is 1.72 for 
EXEM, 2.32 for PL1WM, 2.14 for RDI3WM, 2.17 for RDI4WM, 2.07 for NON 
and order 2.17 for our new scheme DRI1. 

Comparing the computational effort versus precision, in this example the 
schemes DRI1 and RDI4WM perform better than NON and RDI3WM and 
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Fig. 2. Orders of convergence and computational effort per simulation path versus 
precision for SDE (f36j) . 

clearly better than EXEM and the Platen scheme. 

Our last example is a nonlinear SDE with 10 Wiener processes, 



dX{t) = X(t) dt + ^Jx(t) + l - dW.it) + -W*(f) + \ dW 2 {t) 



+yjx(t) + 1 dw 3 (t) + Ux{t) + 1 dw 4 (t) + y X (t) + 1 WW 



+ ^ V x(t) + To dW9{t) + h V x(t) + ^ rfw/lo(t) ' x(0) = L (38; 

Here, we consider the fourth moment, i.e., f(x) = x 4 , and obtain 

4625768169 2998776077847 731453, 80235120932849 251453, . , 

E(f(X(t))) = e 36oooo* + e wow*. (39) 

73570420483600 113706563209000 78178246418000 

We approximate E(f(X(t))) at t = 1 by M = 2 ■ 10 7 simulated trajectories 
with step sizes 2°, . . . , 2~ 3 and obtain the results presented in Table H] and 
Figure El The order of convergence is 1.30 for EXEM, 1.60 for PL1WM, 1.86 
for RDI3WM, 1.91 for RDI4WM, 1.62 for NON and 2.02 for DRI1. If we take 
the computational effort into account, we see that DRI1 performs impressively 
better than all other schemes, which is what we expected for high numbers of 
Wiener processes. 



6 Conclusion 



In the present work, a full classification of the coefficients for a new class of 
efficient explicit SRK methods of order (1, 1) for s — 1 and order (2, 1) for 
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Table 4 

Mean errors, empirical variances and confidence intervals for SDE (|38p . 





h 




*2 


a 


b 




2 U 


-2.793E+01 


7.005E-04 


-2.794E+01 


-2.792E+01 


Jl/VYJl/lVI 


2- 1 


-1.420E+01 


2.521E-03 


-1.421E+01 


-1.418E+01 


2- 2 


-5.658E+00 


7.216E-03 


-5.687E+00 


-5.629E+00 




2- 3 


-1.872E+00 


1.040E-02 


-1.907E+00 


-1.837E+00 




2 -u 


-2.266E+01 


1.183E-03 


-2.268E+01 


-2.265E+01 


r 1j 1 VV 1V1 


2- 1 


-9.218E+00 


1.954E-03 


-9.234E+00 


-9.203E+00 


2- 2 


-2.965E+00 


4.226E-03 


-2.987E+00 


-2.942E+00 




2- 3 


-8.294E-01 


4.294E-03 


-8.519E-01 


-8.070E-01 




2 U 


-1.019E+01 


1.727E-03 


-1.021E+01 


-1.018E+01 


LxlJ LO VV 1V1 


2- 1 


-3.161E+00 


2.324E-03 


-3.177E+00 


-3.144E+00 


2" 2 


-8.582E-01 


4.494E-03 


-8.812E-01 


-8.353E-01 




2-3 


-2.136E-01 


4.373E-03 


-2.363E-01 


-1.910E-01 




2 U 


-9.546E+00 


1.930E-03 


-9.561E+00 


-9.531E+00 


1 V IJlH: VV 1V1 


2- 1 


-2.824E+00 


2.436E-03 


-2.840E+00 


-2.807E+00 




2- 2 


-7.398E-01 


4.557E-03 


-7.629E-01 


-7.167E-01 




2-3 


-1.791E-01 


4.392E-03 


-2.017E-01 


-1.564E-01 




2 U 


5.331E+00 


4.219E-03 


5.309E+00 


5.353E+00 


NON 


2- 1 


1.883E+00 


5.097E-03 


1.858E+00 


1.907E+00 


2- 2 


5.877E-01 


2.975E-03 


5.690E-01 


6.063E-01 




2-3 


1.850E-01 


2.832E-03 


1.668E-01 


2.033E-01 




2° 


-9.465E+00 


1.103E-03 


-9.476E+00 


-9.453E+00 


DRI1 


2- 1 


-2.743E+00 


3.070E-03 


-2.762E+00 


-2.724E+00 


2- 2 


-6.834E-01 


2.531E-03 


-7.006E-01 


-6.662E-01 




2-3 


-1.425E-01 


2.704E-03 


-1.603E-01 


-1.247E-01 




Fig. 3. Orders of convergence and computational effort per simulation path versus 
precision for SDE (f38j) . 

s = 2 stages as well as for order (2, 2) with s = 3 stages is calculated. Based 
on this classification, coefficients for an extension of the deterministic RK32 
scheme to the stochastic case with minimized error constant are given. For 
three examples, this scheme is finally compared with the order two Platen 
and extrapolated Euler scheme, the schemes RD13WM and RDI4WM and 
NON. It turns out that the new developed scheme performs very well and 
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especially much better than all other schemes in the case of a high number of 
Wiener processes. 
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